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Abstract 


We present the basic ideas of the dynamics system of 
the finite-volume General Circulation Model developed at 
NASA Goddard Space Flight Center for climate simula- 
tions and other applications in meteorology. The dynam- 
ics of this model is designed with emphases on conser- 
vative and monotonic transport , where the property of 
Lagrangian conservation is used to maintain the physi- 
cal consistency of the computational fluid for long-term 
simulations. As the model benefits from the noise-free so- 
lutions of monotonic finite-volume transport schemes , the 
property of Lagrangian conservation also partly compen- 


sates the accuracy of transport for the diffusion effects 
due to the treatment of monotonicity . By faithfully main- 
taining the fundamental laws of physics during the com- 
putation , this model is able to achieve sufficient accuracy 
for the global consistency of climate processes. Because 
the computing algorithms are based on local memory , this 
model has the advantage of efficiency in parallel computa- 
tion with distributed memory. Further research is yet de- 
sirable to reduce the diffusion effects of monotonic trans- 
port for better accuracy , and to mitigate the limitation 
due to fast-moving gravity waves for better efficiency. 
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1. Introduction 

The activities of the atmosphere can be simulated with 
computers by programming approximate solutions to the 
physical laws that govern the atmosphere, which is repre- 
sented by a set of grid boxes in the computers. An atmo- 
spheric model is usually divided into two parts: a dynam- 
ical core that solves a set of partial differential equations 
for an idealized atmosphere, and a set of physical param- 
eterizations that provide realistic forcing to the idealized 
atmosphere, such as the precipitation in the atmosphere, 
the friction on Earth’s surface and the radiation from the 
Sun. These two parts are usually treated separately, and 
we will focus on discussing the numerical methods that 
are employed to construct approximate solutions for the 
idealized atmosphere, where there are no sources or sinks 
of physical quantities implied by the physical parameter- 
izations. 

As distributed memory has become a trend of mod- 
ern computer design, the efficiency of the implementa- 
tion of parallel computation with distributed memory 
has become an important issue for designing computing 
algorithms for climate modeling which requires massive 
computation for long-term simulations. While numerical 
methods based on local memory have the advantage of 
parallel efficiency for distributed memory, there are fun- 
damental issues concerning the application of these local 
methods to global climate modeling. 

One of the key concerns for adopting local methods 
for climate modeling is the maintenance of global consis- 
tency among the climate phenomena. To achieve global 
consistency without using global information to constrain 
the numerical solutions, it is necessary for a local method 
to be as faithful as possible to the conservation laws of 
physics. The finite-volume method , which approximates 
the integral of a conservation law over the volume of each 
grid box, is thus considered viable for successful climate 
modeling because of its useful property of local conserva- 
tion for maintaining physical consistency in the numerical 
solutions. Furthermore, the property of Lagrangian con- 
servation, which conserves the physical quantities con- 
tained in a fluid element following the motion, can be 
used to improve the accuracy of the numerical solutions. 
Note that Lagrangian conservation is a special case of 
local conservation, and a physical quantity contained in 
a fluid element is conservative following the motion only 
when there are no sources or sinks during the transport 
process. 

Another general issue of simulating atmospheric mo- 
tions with computers is that spurious noise is often gen- 
erated in the approximate numerical solutions. As spu- 
rious noise contaminates the simulation with unrealistic 


features, various techniques have been developed to miti- 
gate the problem, yet it may be more desirable not to have 
the spurious noise generated in the first place. The spu- 
rious noise can indeed be prevented during the transport 
process by maintaining the monotonicity of the distribu- 
tions of physical quantities in the sense that no new local 
maxima or minima are created by the numerical approxi- 
mation. Unfortunately, monotonic transport schemes are 
often quite diffusive and therefore not necessarily more 
accurate than non-monotonic transport schemes in all 
cases. However, for scalar transport, in the absence of 
sources and sinks, the continuous governing equations 
dictate that the solutions be monotonic and it may be 
a desirable feature that the discrete model maintain this 
property in the solutions. 

Based on the above considerations, we have developed 
a General Circulation Model (GCM) at the Goddard 
Space Flight Center of National Aeronautics and Space 
Administration (NASA), using monotonic finite- volume 
transport schemes with the choice of physical quantities 
that are conservative following the motion in the ide- 
alized smooth atmosphere. We have adopted a set of 
finite-volume schemes for the foundation of transport in 
one dimension, including the piecewise constant scheme 
of Godunov [1], a piecewise linear scheme of van Leer [2] 
and the Piecewise Parabolic Method (PPM) of Colella 
and Woodward [3]. The dynamical core of this model is 
based on the work of Lin and Rood [4, 5, 6, 7] in 1990s, 
and the physical parameterizations are based on those of 
the Community Climate Model (CCM) of the National 
Center for Atmospheric Research (NCAR). In this arti- 
cle, we present the basic ideas of the NASA finite- volume 
dynamical core, concentrating on how to use the prop- 
erty of Lagrangian conservation to improve the accuracy 
of the approximate solutions for the idealized atmosphere. 


2. The transport algorithms 

The status of an idealized atmosphere can be uniquely 
specified with a set of state variables such as the air mass 
density, the air temperature and the wind velocity. Start- 
ing with a set of data of the state variables for the grid 
boxes that represent the idealized atmosphere, the ob- 
jective of a numerical method is to predict the values of 
the state variables for the grid boxes in a short period of 
time, say 5-60 minutes, and a simulation is made with a 
sequence of such small time steps. 

For the purpose of illustration, let us consider the PPM 
transport mechanism of the air in one dimension as shown 
in Fig. 1. Suppose we have the mean density of air mass 
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over the grid intervals at an initial time, and we wish 
to determine the mean density of air mass over the grid 
intervals after a short time of transport. We first notice 
that the air mass contained in a grid interval is simply the 
mean density times the length of the grid interval, and 
the mean density after transport is readily determined if 
we know the amounts of inflow and outflow of the grid 
interval during the transport. We also notice that the 
outflow of a grid interval across a grid interface is exactly 
the inflow of the adjacent grid interval across the same 
interface, they are the same amount of air crossing that 
interface during the transport, i.e., the air mass flux. The 
issue is then down to the determination of the air mass 
fluxes across the grid interfaces during the transport, and 
the accuracy of the approximate solution to the mean 
density after transport depends solely on how the air mass 
fluxes are calculated. 

The PPM transport algorithm starts with the construc- 
tion of a detailed structure of air mass density within each 
grid interval, using the given mean values of air mass 
density for the grid intervals. These detailed structures 
within the grid intervals will be called the subgrid dis- 
tributions and the collection of the mean values will be 
called the overgrid distribution here for convenience of 
discussion. Figure 1 illustrates an example of overgrid 
distribution with black line segments and the correspond- 
ing subgrid distributions with pieces of red curves. Note 
that the subgrid distribution within a grid interval is con- 
structed in such a way that the area under the curve is 
equal to the rectangular area under the line segment. In 
other words, the integral of a subgrid distribution over 
the corresponding grid interval is equal to the air mass 
contained in this grid interval. We also note that each 
subgrid distribution is monotonic within the correspond- 
ing grid interval, and adjacent subgrid distributions are 
monotonic across the grid interfaces. The monotonicity 
of the subgrid distributions is a consequence of the mono- 
tonicity of the given overgrid distribution. 

Now, assume that the velocity of the air flow is known 
everywhere, then we know exactly which air particle 
crosses a grid interface during the transport. Suppose 
the air flow moves from the grid interval A to the grid 
interval B in Fig. 1, then we know how far upstream in 
grid interval A would an air particle reach the interface Ii 
between A and B, and the air contained in this upstream 
range is exactly the flux across interface Ii during the 
transport. In other words, the flux F x across interface Ii 
is equal to the volume integral of the subgrid distribution 
over the upstream range in grid interval A, which corre- 
sponds to the shaded area under the curve in grid interval 
A in Fig. 1. 


Similarly, we can determine the flux F 2 across the in- 
terface I 2 between the grid intervals B and C in Fig. 1. 
Suppose the air flow moves from B to C, then the flux F 2 
is an outflow leaving B during the transport. Given the 
air mass M contained in grid interval B before transport, 
we see that (M - F 2 ) is the amount of air remaining in 
B after transport, i.e., the volume integral of the subgrid 
distribution over the remaining section of B, and it corre- 
sponds to the blank area under the curve in grid interval 
B in Fig. 1. With F\ being the inflow and F 2 the out- 
flow during the transport, the air mass contained in grid 
interval B after transport is simply AT = M + Fi - F 2 , 
and the mean density over B after transport is obtained 
by dividing the resulted air mass AT with the length of 
the grid interval. 

Note that the overgrid distribution of the mean values 
after transport is monotonic if the overgrid distribution of 
the mean values before transport is monotonic, the mono- 
tonicity is maintained by the monotonic construction of 
subgrid distributions. Furthermore, each mean value af- 
ter transport is obtained by combining the neighboring 
volume integrals of the subgrid distributions that con- 
serve the material within the corresponding grid intervals, 
the transported material is thus conserved locally. Since 
these constructed subgrid distributions are not exactly 
the actual ones for the smooth fluid to be simulated, er- 
rors have been introduced during the numerical transport. 
These errors would be greater if the transported material 
were not conservative following the motion in the smooth 
fluid system, for the actual distribution would have been 
changed during the transport-yet the approximate sub- 
grid distributions were constructed with the overgrid dis- 
tribution given at the initial time. 

The errors in one-dimensional transport could be fur- 
ther amplified in multi-dimensional transport when using 
the one-dimensional schemes to estimate the fluxes in 
each dimension separately. Lin and Rood [4] designed a 
multi-dimensional transport algorithm to take advantage 
of the efficiency of dimensional splitting while reducing 
the associated errors by considering the contributions 
from other dimensions during the transport. This multi- 
dimensional algorithm also maintains the relationship 
between air mass and other transported physical quan- 
tities, and thereby promotes the physical consistency of 
the computational fluid system. To further improve both 
accuracy and efficiency in multi-dimensional transport, 
the Lagrangian coordinates [7] is adopted in the vertical 
direction, namely, the altitudes of horizontal material 
surfaces which evolve during the transport. 
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Figure 1: PPM transport mechanism in one dimen- 
sion. Detailed subgrid distributions (red curves) are 
first constructed from the given overgrid distribu- 
tion of mean values (black line segments) for the grid 
intervals A, B and C. The fluxes crossing the grid 
interfaces during the transport are then calculated 
with volume integrals (shaded areas) of the subgrid 
distributions over the upstream ranges implied by 
the given velocity field (green arrows). 

3. The dynamics system 

As illustrated in the previous section, the property of La- 
grangian conservation in the idealized atmosphere can be 
used to improve the accuracy of the approximate solu- 
tions by reducing the errors in the numerical transport 
processes. Not all physical quantities in the continuous 
dynamics system of the atmosphere are conservative fol- 
lowing the motion. For instance, the momentum (mass 
times velocity) carried in an air parcel is not conservative 
following the motion because of the general existence of 
external forcing to the air parcel, such as the pressure- 
gradient force due to the pressure differences surrounding 
the air parcel. If the momentum were chosen to be con- 
served in the discrete numerical system through explicit 
transport with the wind velocity, the errors in inertial 
transport would be amplified by the interaction with the 
pressure-gradient force, and the conservation of momen- 
tum in the discrete system would not be localized as well 
as in the inertial case. 

There are three physical quantities in the continuous 
dynamics system of the atmosphere that are conservative 
following the motion, they are the air mass , the potential 
temperature and the absolute vorticity. Potential tem- 
perature is a measure of molecular motions with respect 
to a reference thermodynamic state of the atmosphere. 
Absolute vorticity is a measure of the rotation of a fluid 
element relative to an inertial frame such as the space 
containing the Earth. These three quantities have been 


chosen to be conserved in the NASA finite-volume dy- 
namical core, to promote the accuracy in transport and 
the localization of conservation. 

Because they are directly related to primary variables, 
the properties of local conservation and monotonicity of 
air mass and potential temperature can be achieved di- 
rectly by finite-volume transport. Absolute vorticity is, 
however, a quantity derived from the primary variable ve- 
locity. To achieve the local conservation and monotonic- 
ity of absolute vorticity with velocity being a primary 
variable, it is essential to apply the Circulation Theorem 
that the integral of vorticity over an area is equal to the 
integral of velocity along the boundary enclosing the area. 
By using the so-called D-grid configuration as shown in 
Fig. 2a, the primary variable velocity can be predicted 
directly along the boundary of the area while conserving 
absolute vorticity. Lin and Rood [5] demonstrated the ad- 
vantage of such an approach in two dimensions with an 
excellent simulation of the Rossby-Haurwitz waves which 
has been recognized as a standard test for modeling at- 
mospheric dynamics [8]. 

To extend the two-dimensional dynamics to three 
dimensions, we have adopted Lagrangian coordinates 
[7] in the vertical direction. Note that the vertical 
Lagrangian coordinates are the altitudes of horizontal 
material surfaces, and we may assume the hydrostatic 
equilibrium that the weight of the air confined between 
two material surfaces are balanced by the difference of the 
pressure at the material surfaces. The three-dimensional 
atmosphere is thus decoupled into a sequence of inertial 
horizontal layers, and the dynamics in each layer can be 
treated in a manner very similar to the two-dimensional 
dynamics, except the calculation of pressure-gradient 
force. The pressure- gradient force in a Lagrangian 
layer [6] is calculated with the Green Theorem that the 
integral of pressure gradient over the volume contained 
in a grid box is equal to the integral of pressure over the 
surface surrounding the grid box as shown in Fig. 2b. In 
other words, the net external force acting on a fluid ele- 
ment is obtained by integrating the surrounding pressure. 

(a) (b) 
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Figure 2: (a) D-grid configuration for applying the 
Circulation Theorem in two dimensions. The abso- 
lute vorticity (ft) is a mean value over the grid box, 
and the wind components (u, v) are predicted along 
the boundary, (b) Schematic for applying the Green 
Theorem with Lagrangian coordinates in vertical di- 
rection. The net external force acting on a fluid 
element is obtained by integrating the surrounding 
pressure (blue arrows). 

The NASA finite-volume dynamical core is in the end 
a collection of finite volumes where dynamical processes 
are calculated as integrals of simple locally continuous 
functions. Many of the inherent errors of finite-difference 
and spectral methods are reduced, and there is a high 
degree of consistency of primary variables and dynamical 
quantities. Figure 3 demonstrates the climatology of 
the NASA finite-volume GCM with the 500-hPa eddy 
height, i.e., the height at 500-hPa pressure level relative 
to its own zonal mean, for December, January and 
February. The simulation is produced with a 15-year 
run at the resolution of 2.5 degrees in longitude and 
2.0 degrees in latitude, and the vertical dimension of 
the atmosphere is covered with 55 Lagrangian layers 
up to the pressure level of one Pascal. The size of 
time step is 7.5 minutes for the dynamical core, and 
it is 30 minutes for the physical par ameterizat ions. 
The PPM is adopted for the foundation of transport 
in one dimension. Note that the 500-hPa eddy height 
represents the wave patterns and magnitudes in middle 
troposphere, and the simulated climatology is quite 
realistic according to the analysis from the European 
Center for Medium-range Weather Forecast (ECMWF). 
Additional results from this model are available at the 
web site <http://dao.gsfc.nasa.gov/NASanCAR> . 


4 • Future development 

Our experience with the NASA finite-volume GCM indi- 
cates that the diffusion effects of the adopted monotonic 
transport schemes are often too strong for simulating 
small-scale phenomena in the atmosphere. We have de- 
veloped new monotonic finite-volume transport schemes 
to address this issue, some of these new transport schemes 
have been incorporated in the model as options, some are 
yet to be published and implemented. 

While the property of Lagrangian conservation can be 
used to enhance the physical consistency in the numeri- 
cal solutions, the accuracy of transport depends largely 
on how the velocity is predicted in the dynamics system. 
Although Lin and Rood [5,6] provide an accurate way 
to predict the velocity while taking advantage of the La- 
grangian conservation of absolute vorticity, the efficiency 
of their methodology is subject to the restriction on the 
size of time steps for the fast-moving gravity waves. Fur- 
ther research is desirable to enlarge the time step and 
hence improve the overall efficiency. 

Furthermore, due to the adoption of the traditional 
longitude-latitude grid system, global memory is still 
needed in polar areas where the meridians converge, and 
this affects the efficiency of the implementation of paral- 
lel computation with distributed memory. A new design 
of the dynamical core with quasi-uniform grids is under- 
going development, in order to refrain from using global 
memory, and thereby promote the parallel efficiency for 
distributed memory. It is yet to be proved that this ap- 
proach would gain enough in efficiency while achieving 
sufficient accuracy for practical applications. 

For more information about the NASA finite-volume 
GCM, the readers are referred to the documentation 
on the next-generation model provided at the web site 
<http://dao.gsfc.nasa.gov/pages/atbd.html>. 


(this figure is provided in a separated file) 


Figure 3: Climatology simulated with the NASA 
finite-volume GCM (upper panel) and the corre- 
sponding analysis from ECMWF (lower panel), 
demonstrated with an average of 500-hPa eddy 
height (m) over December, January and February. 
The simulation is produced with a 15- year run at 
the resolution of 2.5x2.0 degrees with 55 layers in 
vertical direction. 
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The dynamics of this model is designed with emphases on conservative and 
monotonic transport, where the property of Lagrangian conservation is used 
to maintain the physical consistency of the computational fluid for long-term 
simulations. As the model benefits from the noise-free solutions of monotonic 
finite-volume transport schemes, the property of Lagrangian conservation also 
partly compensates the accuracy of transport for the diffusion effects due to 
the treatment of monotonicity. 

By faithfully maintaining the fundamental laws of physics during the 
computation, this model is able to achieve sufficient accuracy for the 
global consistency of climate processes. 

Because the computing algorithms are based on local memory, this model has 
the advantage of efficiency in parallel computation with distributed memory. 

Further research is yet desirable to reduce the diffusion effects of 
monotonic transport for better accuracy, and to mitigate the limitation 
due to fast-moving gravity waves for better efficiency. 


